Convergence and stability in numerical relativity 
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It is often the case in numerical relativity that schemes that are known to be convergent for well 
posed systems are used in evolutions of weakly hyperbolic (WH) formulations of Einstein's equations. 
Here we explicitly show that with several of the discretizations that have been used through out the 
years, this procedure leads to non-convergent schemes. That is, arbitrarily small initial errors are 
amplified without bound when resolution is increased, independently of the amount of numerical 
dissipation introduced. The lack of convergence introduced by this instability can be particularly 
subtle, in the sense that it can be missed by several convergence tests, especially in 3+1 dimensional 
codes. We propose tests and methods to analyze convergence that may help detect these situations. 
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Convergence is a central element of any numerical sim- 
ulation. It refers to the property that if one refines the 
simulation by adding more points to the grid, numerical 
errors should diminish. In the limit of zero spacing, they 
should go to zero and one should get the exact solution. 
Having a convergent code is a key element for numerical 
simulations to have predictive power: although one will 
always be limited in practice to a finite number of points 
in the grid, one can extrapolate the results for more and 
more refined simulations and have a very good approx- 
imation to the true (continuum) results. Codes that do 
not converge produce answers that, even if they remain 
finite (at least for a while), have little predictive power: 
there generically is no way to know if the solutions found 
approximate the desired continuum solution. 

In this paper we want to emphasize that discretization 
schemes that yield convergent code for strongly hyper- 
bolic (SH) systems of equations do not necessarily do 
so for weakly hyperbolic (WH) or completely ill posed 
(CIP) systems. The relevance of this observation is 
that several formulations of the Einstein evolution equa- 
tions commonly used in numerical relativity, including, 
for example, the Arnowitt-Deser-Misner (ADM) [ 1 ] and 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) [2] for- 
mulations with fixed lapse and shift, are not SH. 

We concretely prove that a simple system of linear WH 
equations with constant coefficients, when discretized 
with either the iterated Crank-Nicholson (ICN) with 
fixed number of iterations or second order Runge-Kutta 
(2RK) methods, leads to unconditionally unstable codes, 
even if numerical dissipation is explicitly added. This 
simple system of equations is related to the equations 
one encounters in the WH formulations of general relativ- 
ity. It therefore strongly suggests that these formulations 
should produce code that is not convergent. The lack of 
convergence is of a particularly pernicious nature, since 
in the WH case it might grow slower than the "usual" 
von-Neumann numerical instability. It is such that if 
one tests the code with stationary solutions (for instance 
simulating a single black hole) and performs convergence 
tests, one could easily be confused into believing one has 
convergent code. 

Consider a linear system of m partial differential equa- 



tions in m variables of the form 

d t v = Ad x v , (1) 

where A is a constant m X m matrix and v = 
(v\, . . . , v m ) T . The system is SH if A has real eigenvalues 
and is diagonalizable, WH if it has real eigenvalues but is 
not diagonalizable, and CIP if it has complex eigenvalues. 

The iterated Crank-Nicholson (ICN) method is an ex- 
plicit approximation to the (implicit) Crank-Nicholson 
method first proposed by Choptuik. In this proposal 
the number of iterations is not fixed but, instead, might 
depend on the spatial gridpoint and time step A 
different notion of ICN is to a priori fix the number of 
iterations |§. Teukolsky || showed that for the advec- 
tive equation, dtv = d x v (which is SH), such scheme was 
stable if one used 2, 3, 6, 7, 10, 11, . . . iterations (one iter- 
ation corresponds to the 2RK method). This paper seeks 
to extend this proof to more general systems of equations, 
including explicit artificial viscosity. 

A difference scheme is said to be stable with respect 
to a given norm if there exist positive constants Aa;o 
and At and f(t) such that \\u n \\ < f(t)\\u°\\, for 
< Ax < Ax , < At < Ato, and t = nAt where 



is the numerical solution at time t and 



the 



tial data. That is, the growth of the solution is bounded 
by some function of time fit) that is independent of the 
initial data and resolution. Some schemes are only condi- 
tionally stable. This means that a relationship between 
At and Ax is needed for stability. Through Lax's the- 
orem, numerical stability in this sense is equivalent to 
convergence, provided the scheme is consistent. 

Working in Fourier space and writing explicitly the so- 
lution at time t in terms of the amplification matrix Q, 
u n = Q n u° yields a necessary condition for stability: the 
von-Neumann condition (VNC). It states that, for the 
cases considered in this paper, the spectral radius p{Q) 
of Q (the maximum eigenvalue in norm) is not greater 
than one. If this condition is not satisfied, a numeri- 
cal instability that grows like exp(n/i), with /i > a 
constant, is present. A well known example of this in- 
stability is the forward in time, centered in space scheme 
for the advective equation. It is sometimes thought that 
the VNC is not only necessary but also sufficient for sta- 
bility. This often leads to (wrongly) concluding that a 
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discretization for a WH system is stable based on an dis- 
crete analysis for the advection or wave equation. The 
VNC is sufficient only if Q can be diagonalized by a simi- 
larity transformation whose norm, and that of its inverse, 
has an upper bound that is independent of resolution [Q . 
An observation that will be important for the stability 
analysis below is that, for the cases we are discussing, Q 
is diagonalizable if and only if A is. Therefore, in the 
WH case the VNC will not be sufficient for stability. 

On the other hand, if the system is SH A is diago- 
nalizable (with real eigenvalues) and so is Q. Working 
in the diagonal basis we then have an uncoupled sys- 
tem of to equations for the to grid functions. Repeating 
Teukolsky's analysis for the ICN scheme, this system is 
stable if p(A)X < 2 (A = At /Ax is the Courant factor) 
and the number of iterations is p = 2, 3, 6, 7, 10, 11, . . ., 
if no explicit dissipation is added. In evolutions of SH 
equations, the addition of dissipation may help stabilize 
schemes that are otherwise unstable. For example, the 
VNC for the 2RK case with third order dissipation is 

(p{A)X) i < 8\a < 1 , (2) 

where a is the dissipation parameter. Since Q is diagonal, 
this condition is necessary and sufficient for stability. 

Below we will explicitly prove that for CIP or WH cases 
the ICN and 2RK methods are unstable, independently 
of the amount of dissipation, number of iterations and 
Courant factor. Although our results are only proven for 
a system of linear equations with constant coefficients, 
they are significant in several ways for numerical relativ- 
ity. First of all, if one considers situations that are small 
perturbations off Minkowski spacetime, then a system of 
equations identical to the one we considered does appear 
in the linearized Einstein equations written in first order 
form in space and time if the system of evolution equa- 
tions used is WH or CIP. Since far away from binary black 
holes the spacetime is approximately flat, this situation 
does arise in realistic simulations. Could non-linearities 
change our conclusions? It is possible, but unlikely. Usu- 
ally when instabilities are established for linear systems 
the addition of non-linear terms only makes instabilities 
worse 0. 

The type of lack of convergence differs in the CIP and 
WH case. In the CIP case the VNC is easily violated, 
codes crash very fast and non-convergence is easy to spot. 
In the WH case the VNC may be satisfied and crashes 
may occur very late (if they occur at all). If one per- 
forms routine convergence analyses in the region before 
the crash one might appear to see convergence, especially 
if the initial data has few (and low) frequency compo- 
nents, unless one increases enough the resolution or runs 
for long enough time. 

We now show that these different types of lack of con- 
vergence do appear in numerical relativity simulations, 
even in very simple ID ones with periodic boundary con- 
ditions. We have performed non-linear evolutions of a 
plane gauge wave spacetime (that is, the spacetime is 



flat, though in a time-dependent slicing), 

ds 2 = e Ash ^ t+x \-dt 2 + dx 2 ) + dy 2 + dz 2 . (3) 

Figure [l] shows results of an evolution using 2RK with 
a = 1/2,0, -1/2 in the Einstein-Christoffel (EC) system 
H of evolution equations. The latter is a symmetric- 
hyperbolic reformulation of Einstein's evolution equa- 
tions that includes a parameter a that densitizes the 
lapse. In the EC formulation, a = 1/2, but by tuning 
this parameter we can make the system WH (a — 0) or 
CIP (a = -1/2). The whole 30 equations of the EC 
system are evolved, but all quantities are assumed to de- 
pend only on t and x s [-— 7r,7r]. The Courant factor is 
set to 1/2, the dissipation parameter to a = 0.02, and 
the evolution is followed for 1, 000 crossing times. In the 
SH case the code is convergent for all resolutions tested, 
see j(| for a detailed discussion. On the other hand, in 
the same code with a = —1/2 a lack of convergence be- 
comes apparent immediately (before one crossing time), 
the errors become bigger and the code crashes earlier as 
resolution is increased. In the WH case the code is also 
not convergent but the effect is less noticeable, since with 
the chosen values of dissipation and Courant factor the 
VNC, Eq.(|), is satisfied. 

If one performed only a few runs, with, say 120 and 
240 gridpoints as is typical in 3D convergence studies, 
one would have to wait for around 150 crossing times in 
order for the lack of convergence in this WH example 
to become obvious. To put these numbers in context, 
suppose one had a similar situation in a 3D black hole 
evolution. Suppose the singularity is excised, with the 
inner boundary at r = M and the outer boundary at 
20M. In this case 120 and 240 gridpoints correspond to 
grid spacings of, approximately, M/5 and M/10, respec- 
tively. If one had to evolve up to 150 crossing times in 
order to notice the lack of convergence, that would cor- 
respond tot k 3000 M, which is more than what present 
3D evolutions typically last. One could therefore be mis- 
led to think that the code is convergent. Repeating the 
same runs with an initial data with more frequencies (say, 
a non-stationary black hole) would make the instability 
manifest in a shorter timescale. By making a spatial 
Fourier decomposition (Figure |^) of the numerical so- 
lution, we have found that in the WH and CIP cases 
there are always non-zero frequency modes growing ex- 
ponentially from the very beginning, though sometimes 
starting at truncation error. By performing such decom- 
position one can detect that the code is not convergent 
much before this becomes obvious in the overall errors. 

We have done simulations with different spectral distri- 
butions on the initial data, different number of iterations 
in the ICN method, different Courant factors, and dif- 
ferent values of dissipation. The time and resolution at 
which the numerical instability becomes obvious depends 
on all these factors, but lack of convergence is always 
present in the CIP and WH cases, while the SH runs do 
converge. Too much dissipation violates the VNC, result- 
ing in a more severe instability and this is immediately 
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FIG. 1: L2 norm of the errors for the metric in the SH, WH 
and CIP cases. 



seen in numerical experiments, see Q for details. 

We have also found similar results performing simula- 
tions with the same initial data but with Kidder-Scheel- 
Teukolsky's many parameter family of formulations of 
Einstein's equations , leaving all parameters fixed and 
changing only one of them at a time, achieving different 
levels of hyperbolicity. In particular, we have found lack 
of convergence in the ADM equations as well (which is a 
subsytem of this family). 

The results presented in this paper do not represent a 
special feature of the ICN method but instead, only re- 
flect the fact that the definition of numerical stability is 
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FIG. 2: The lack of convergence in the particular WH case 
considered is triggered by non-zero frequency modes that are 
suppressed in the initial data, but grow exponentially. Notice 
that after the lack of convergence becomes obvious (at t = 
2000 for this resolution) no further exponential growth takes 
place and the code does not crash, in the example we are 
considering. 



just a discrete version of well posedness. Therefore, dif- 
ference schemes approximating ill posed problems can be 
expected to be non-convergent. In the absence of bound- 
aries this is the case for CIP or, generically, WH formla- 
tions. If there are boundaries, strong hyperbolicity is not 
enough and extra care has to be taken in order to guaran- 
tee well posedness; wrong boundary conditions not only 
lead to inconsistencies but also to lack of convergence 
(see, e.g., [§). 

Although it is not possible to prove in a definitive way 
that a code is convergent only by numerical experiments, 
the previous examples suggest some of the pathologies 
that one should look for. Namely, Fourier modes that 
are not convergent but that are very small for a while, 
since they are initially suppressed. Also, notice from the 
WH simulations above shown that a code does not need 
to crash nor exhibit violent growth in order not to con- 
verge. The main lesson learned from this paper is that 
one should exercise significant care in numerical simula- 
tions before empirically concluding that the simulation 
is convergent, especially if the formulation of Einstein's 
equations used is not SH or its level of hyperbolicity is 
unknown. 

We wish to thank D. Arnold, M. Choptuik, P. Laguna, 
L. Leaner, M. Miller, R. Price, O. Rcula, and S. Teukol- 
sky for comments. This work was supported in part by 
NSF grant PHY9800973, the Horace Hearne Jr. Insti- 
tute for Theoretical Physics, the Swiss National Science 
Foundation, and Fundacion Antorchas. 

Proof of non-convergence for CIP and WH cases: A 
usual way of discretizing the right hand side of Eq. ([j]) is 
using centered differences plus third order explicit dis- 
sipation. By this one means solving dtv = Cv/At, 
with C = AXSq/2 — (iA<5 4 , where I is the identity ma- 



4 



trix, S v k = v k+1 - v k -i, and 8 A v k = v k+2 - 4u fe+ i + 
6v k — 4:V k _i + v k -2- This is equivalent to discretizing 
d t v — Ad x v — Ia(Ax) 3 d x v using centered differences (us- 
ing first order dissipation in our proof, i.e. discretizing 
dtv = Ad x v + dAxd^v, gives similar results). 

The ICN method with p iterations for the partial dif- 
ferential equation (til) can be written as 



with e arbitrary small. In Fourier space we have u° (£) = 
e(0, 1) T for all £. The solution at time t = nAt is, 
then, u" = e(a n ~ 1 nb, a n ) T . We will show that its norm 
grows without bound when the number of gridpoints is 
increased (while keeping the time and Courant factor 
fixed). We first notice that (712 = 2n — 2), 
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Expressions |a| 2 and |6| 2 are analytic functions of £ with 
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The index n corresponds Taylor expansion \b\ 2 = A 2 [£ 2 + 0(£ 4 )] and 




to the time step and k to the spatial mesh point. Since 
we are considering an initial value problem on all of IR, 
— 00 < k < +00 (our proof can be easily modified for pe- 
riodic boundary conditions, and the results are the same). 

In order to analyze stability we work in the basis in 
which A takes its canonical form. That is, one multiplies 
both sides of Eq.(g) by T, where TAT^ 1 = J has the 
canonical Jordan form, and analyzes the equation 



(5) 



where L = TCT^ 1 = J\S /2 - IZXS 4 . Any con- 
clusions regarding stability will hold as well for the 
original variable v. Working in Fourier space, u k = 
£ n exp Q = 1 + ZTijtlti/V and L = 

JM sin(0 - 16ffAsin 4 (£/2). 

If the system is CIP, J has at least one complex eigen- 
value c, say Jn = c. In this case, Ln = cAisin(£) — 
16(tA sin 4 (£/2), and the eigenvalue Qn has norm 1 — 
XIm(c)£ + 0(£ 2 ). Therefore, the VNC is violated for suf- 
ficiently small £ with the appropriate sign and the scheme 
is unconditionally unstable, as expected. 

If the system is WH then J has at least one Jordan 
block and so does Q. In this basis all the Jordan blocks 
are uncoupled and we can consider one at a time. We 
assume there is one block of dimension 2x2 (the proof 
actually holds for higher dimensionality as well). The 
canonical form for such a block, the resulting amplifica- 
tion matrix, and its n-th power are 
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Suppose the VNC \a\ < 1 is satisfied (this rules out the 
p = case), we will show that the scheme is still unstable. 
Consider as initial data a small perturbation of the form 
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(0, 0) T otherwise (6) 



= 1 - 

[2a\ + d(cA) 4 ]£ 4 + 0(<f ), with d = -1 for p = 1, d = 1 
for p — 2, and d — otherwise. From these expansions 
one can see that, for all p, there are positive constants a 
and p such that 



H£)| 2 >i-< 4 >o 



H0\ 2 >p 2 f >o, 



(8) 



for small enough f, say |£| < e\. For n > 1 and small 
enough say |£| < e 2 , 



|a(0| 2 " > l-na£ 4 > , 



(9) 



where the last inequality holds for |£| < (no;) -1 / 4 , 
provided that n is large enough so that (net) -1 / 4 < 
min(ei, £2). Using bounds (|[||) and the VNC in inequal- 
ity (0) gives 



,n||2 
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which diverges for n — > 00. This shows that in the WH 
case the ICN and 2RK schemes are unstable, independent 
of the amount of dissipation. In fact, by examining |a(£ = 
7r) I , one can see that a necessary condition for the VNC 
(for any number of iterations) is < <tA < 1/8. Adding 
too much dissipation violates the VNC and the instability 
becomes worse. 

Perturbations like that of (g) are to be expected in a 
numerical simulation due to truncation or roundoff er- 
rors. For high enough resolution, such a perturbation 
will be amplified without bound, spoiling any conver- 
gence. One expects that in the non-constant coefficient 
or in the non-linear case the rate of growth of the insta- 
bility with number of gridpoints will be even faster (see 
example in page 216 of Qj), as in figure |[ 
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